function [out_pic,score] = partition(edge_map, penalty_vec, d)
% partition finds the highest score partition int the edge map using
% dynamic programing.
    
    % intializing input sizes
    if nargin < 2 
        penalty_vec = [0,0,0];
    end;
    if nargin < 3
        d = 0;
    end;

    % disabeling frequency when d = 0
    if  d == 0
            freq_weight = 0;
    else
        freq_weight = 0.3/d;
    end
    
    [N,M] = size(edge_map);
    [~,range] = size(penalty_vec);
    range = (range-1)/2;
    penalty_vec = penalty_vec - max(penalty_vec); % making all the penalties negative (no bonuses) to not afect the top score
    
  %range should be odd, because we are checking in a symmetric manner.
    if mod(range,1) ~= 0
        error('partition:ilegal_args','penalties vector size must be odd ( usualy simetric )');
    end
        
    % creating DP tabels
    dp_table = -Inf(N+2*range,M);
    dp_table(1+range:(N+range),:) = zeros(N,M);
    dp_table(1+range:(N+range),1) = edge_map(:,1);
    ptr_table = zeros(N,M);
    length_table = ones(N+2*range,M);
    phase_table = zeros(N+2*range,M) +d; % contations values from 1 to 2d, depends where are we. 

    % 0 means we don't have preference to any direction
    
    % iterating over the columns
    for j = 2:M
        % foreach line in the column j
        for i = (1:(N))+range 
            % checking foreach cell which is the best option in range
            scores_vec =dp_table(i-range:i+range,j-1);
            phase_vec = phase_table(i-range:i+range,j-1);
            
            % the length of the current track, in case we start with a
            % non vertical location.
            length_vec = length_table(i-range:i+range,j-1);
            
            % we choose the path with the maximum average value, since the
            % paths length may differ, we use length_vec.
            options_vec = (scores_vec + penalty_vec') ;

			% adding the phase factor
            % going down
            options_vec(1:range) = options_vec(1:range) + (phase_vec(1:range) - d) * freq_weight;
            %going up
            options_vec(range+2:2*range+1) = options_vec(range+2:2*range+1) + (d - phase_vec(range+2:2*range+1)) * freq_weight;
            options_vec(range+1) = options_vec(range+1) +0.3* max((phase_vec(1+range) - d) * freq_weight,(d - phase_vec(range+1)) * freq_weight);
            % maybe d-1 here? (line above)
            options_vec = options_vec./ length_vec;
            [val,index] = max(options_vec);
            
            % prefer straight lines instead of sharp transitions when score
            % is equal
            if options_vec(index) == options_vec(range+1)
                index = range+1; 
            elseif options_vec(index) == options_vec(range)
                index = range;
            elseif options_vec(index) == options_vec(range+2)
                index = range+2;
            end
            
                        % update phase table
             if (index <= range)
                 if (phase_table(i+ index - range -1,j-1) >= d)
                     phase_table(i,j) = mod(phase_vec(index) +1 , 2*d);
                 else
                     phase_table(i,j) = d+1;
                 end %if
             elseif (index > range + 1)
                if (phase_table(i+ index - range -1,j-1) < d)
                     phase_table(i,j) = mod(phase_vec(index) + 1, 2*d);
                else
                   phase_table(i,j) = 1;
                end %if
             else %index = range+1
                 phase_table(i,j) = mod(phase_vec(index) + 1, 2*d);
             end %big if
            
            % eliminate transfer from ilegal
            %dp_table(i,j) = edge_map(i-range,j)+scores_vec(index)+penalty_vec(index) ;%+ phase_table(i,j);
            dp_table(i,j)  = val*length_vec(index) + edge_map(i-range,j);
            length_table(i,j) = length_vec(index)+1;
            
            % keeping a reference to the cell the value was taken from
            ptr_table(i-range,j) = (length_table(i,j)>1)*((i- range)+(index-(range+1)));


             
        end
 
    end
    
    % normalizing the length
    dp_table = dp_table./length_table;
    
    % finding the best end position
    [val_right,pos_right] = max(dp_table(:,M));
    
    
    val_top = dp_table(1+range,1:M);
    [val_up,pos_up] = max(val_top);
    
    val_bottom = dp_table(N+range,1:M);
    [val_down,pos_down] = max(val_bottom);
      
    if (val_right >= val_up) && (val_right >= val_down)
             % val_right is max
            pos = pos_right;
            j = M;
            score = val_right;
    elseif (val_up >= val_down)
            % val_up is max
            pos = 1;
            j = pos_up;
            score = val_up;
    else
             % val_down is max
            pos = N;
            j = pos_down;     
            score = val_down;
    end
       
    % marking the path ( finishes at one of the edges or at -N*M )
    out_pic = zeros(N,M);

    while j>0 
        if pos == 0
            break
        end
        out_pic(pos,j) = 1;
        pos = ptr_table(pos,j);
        j = j-1;
    end
   
end
